remove(list=ls())

set.seed(31489)

library(dplyr)
library(ggplot2)
library(tidyr)
library(tidyverse)
library(tibble)
library(haven)
library(Hmisc)
library(AER)
library(doBy)
library(foreign)
library(countrycode)

library(ggeffects)
library(lme4)

library(lmtest)
library(AER)
library(stargazer)
library(fixest)
library(sjPlot)
library(gmodels)
library(MuMIn)

#Set working directory to source file
setwd("C:/APFP_replication")


load("APFP_dat.Rdata")


#Note that descriptives for Table 1 are generated following the regressions for Table 2 below 

#############################
#Models for Table 2##########
#############################

lmer_fat<-lmer(ch_troops ~  lag_troop_deaths + Contiguous + lag_lnpivotalaid + lag_n_contributors_troops + lag_mission_troops + lag_troops_mean + lag_dem + lag_lngdp + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_troops_mean>0, lag_mission_troops>0))
nobs(lmer_fat)
summary(lmer_fat)


lmer_contig<-lmer(ch_troops ~  lag_troop_deaths*Contiguous + lag_lnpivotalaid + lag_n_contributors_troops + lag_mission_troops + lag_troops_mean + lag_dem + lag_lngdp + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_troops_mean>0, lag_mission_troops>0))
nobs(lmer_contig)
summary(lmer_contig)



lmer_aid<-lmer(ch_troops ~  lag_troop_deaths*lag_lnpivotalaid + Contiguous + lag_n_contributors_troops + lag_mission_troops + lag_troops_mean + lag_dem + lag_lngdp + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_troops_mean>0, lag_mission_troops>0))
nobs(lmer_aid)
summary(lmer_aid)


lmer_full<-lmer(ch_troops ~  lag_troop_deaths*lag_lnpivotalaid*Contiguous + lag_n_contributors_troops + lag_mission_troops + lag_troops_mean + lag_dem + lag_lngdp + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_troops_mean>0, lag_mission_troops>0))
nobs(lmer_full)
summary(lmer_full)

r.squaredGLMM(lmer_fat)
r.squaredGLMM(lmer_contig)
r.squaredGLMM(lmer_aid)
r.squaredGLMM(lmer_full)

#Table 2 output (note: standard deviations for random effects and r^2 manually inputted in ms)
stargazer(lmer_fat, lmer_contig, lmer_aid, lmer_full, title="Results from Multi-level Linear Regression", out="table1_replication_6.22pm.doc",
          align=T, digits=2, intercept.bottom=T, model.names=T, dep.var.labels="Change in Troops", covariate.labels=c("Troop Deaths", "Contiguous", "ln(Pivotal Aid)",
                                                                                                                       "No. Contributors", "Contributor Troops",  "Mission Troops", "Democracy", "ln(GDP)", "ln(Battle Deaths)", "Troop Deaths*Contiguous", "ln(Pivotal Aid)*Contiguous",
                                                                                                                       "Troop Deaths*ln(Pivotal Aid)*Contiguous", "Troop Deaths*ln(Pivotal Aid)"), star.cutoffs=c(0.05, 0.01, 0.001), no.space=T, style="io")


########################
#Table 1 Descriptives
########################

stargazer(lmer_full@frame, digits=2, summary.stat=c("mean", "median", "sd", "min", "max"), style="io")




#Generate Figures 1:3

lmer_dat<-model.frame(lmer_full)


contdf<-ggpredict(lmer_full, terms=c("lag_troop_deaths", "Contiguous"))

aiddf<-ggpredict(lmer_full, terms=c("lag_troop_deaths", "lag_lnpivotalaid[1.69, 19.85]"))

fatdf<-ggpredict(lmer_full, terms=c("lag_troop_deaths", "lag_lnpivotalaid[1.69, 19.85]", "Contiguous"))




#Figure 1: Contiguous plot

set_theme(
  title.size=.8,
  axis.title.size=.8,
  axis.textsize=.8,
  legend.size=.6,
  legend.title.size=.7
)


png('fig1_II.png', units="in", width=7, height=3.5, res=300)


plot(contdf, col="bw") + labs(x="Troop Deaths(t-1)", y="Change in Troops") + ggtitle("") + geom_rug(data=lmer_dat, aes(x=lag_troop_deaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()


##Figure 2: Pivotal Aid Plot


png('fig2_II.png', units="in", width=7, height=3.5, res=300)

plot(aiddf, col="bw") + geom_line(aes(linetype=group)) + labs(x="Troop Deaths(t-1)", y="Change in Troops", colour="ln(Pivotal Aid(t-1))") + ggtitle("") + geom_rug(data=lmer_dat, aes(x=lag_troop_deaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()


##Figure 3: Three-way interaction plot
png('fig3_II.png', units="in", width=7, height=3.5, res=300)

plot(fatdf, col="bw") + geom_line(aes(linetype=group)) + facet_wrap(~facet) + labs(x="Troop Deaths(t-1)", y="Change in Troops", colour="ln(Pivotal Aid(t-1))") + ggtitle("") + geom_rug(data=lmer_dat, aes(x=lag_troop_deaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()




####################
#Table 3
#Conflict deaths
###################


lmer_bdeaths<-lmer(ch_troops ~  lag_ln_bdbest*lag_lnpivotalaid*Contiguous+ lag_n_contributors_troops + lag_mission_troops + lag_troops_mean + lag_dem + lag_lngdp + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_troops_mean>0, lag_mission_troops>0))
nobs(lmer_bdeaths)
summary(lmer_bdeaths)

lmer_chbdeaths<-lmer(ch_troops ~  lag_chlnbdbest*lag_lnpivotalaid*Contiguous + lag_n_contributors_troops + lag_mission_troops + lag_troops_mean + lag_dem + lag_lngdp + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_troops_mean>0, lag_mission_troops>0))
nobs(lmer_chbdeaths)
summary(lmer_chbdeaths)

#Additional marginal plots to visualize interactions 

# bdeaths_dat<-model.frame(lmer_bdeaths)
# cont_bdeaths_dat<-ggpredict(lmer_bdeaths, terms=c("lag_ln_bdbest", "Contiguous"))
# 
# plot(cont_bdeaths_dat, col="bw") + labs(x="ln(Battle Deaths(t-1))", y="Change in Troops") + ggtitle("") + geom_rug(data=bdeaths_dat, aes(x=lag_ln_bdbest, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)
# 
# 
# chbdeaths_dat<-model.frame(lmer_chbdeaths)
# cont_chbdeaths_dat<-ggpredict(lmer_chbdeaths, terms=c("lag_chlnbdbest", "Contiguous"))
# 
# plot(cont_chbdeaths_dat, col="bw") + labs(x="Ch. ln(Battle Deaths(t-1))", y="Change in Troops") + ggtitle("") + geom_rug(data=chbdeaths_dat, aes(x=lag_chlnbdbest, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)
# 
# 
# aid_chbdeaths_dat<-ggpredict(lmer_chbdeaths, terms=c("lag_chlnbdbest", "lag_lnpivotalaid"))
# 
# plot(aid_chbdeaths_dat, col="bw") + labs(x="Ch. ln(Battle Deaths(t-1))", y="Change in Troops") + ggtitle("") + geom_rug(data=chbdeaths_dat, aes(x=lag_chlnbdbest, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)
# 
# 
# aidcont_chbdeaths_dat<-ggpredict(lmer_chbdeaths, terms=c("lag_chlnbdbest", "lag_lnpivotalaid[0, 19.85]", "Contiguous"))
# 
# plot(aidcont_chbdeaths_dat, col="bw") + geom_line(aes(linetype=group)) + facet_wrap(~facet) + labs(x="Ch. ln(Battle Deaths(t-1))", y="Change in Troops", colour="ln(Pivotal Aid(t-1))") + ggtitle("") + geom_rug(data=chbdeaths_dat, aes(x=lag_chlnbdbest, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)



#Table 3 output
#Note variable labels, standard deviations for random effects, and r^2 manually inputted in ms manually 

stargazer(lmer_bdeaths, lmer_chbdeaths, title="Results from Multi-level Linear Regression", 
          align=T, digits=2, intercept.bottom=T, dep.var.labels="Change in Troops", star.cutoffs=c(0.05, 0.01, 0.001), no.space=T, style="io")

r.squaredGLMM(lmer_bdeaths)
r.squaredGLMM(lmer_chbdeaths)



###################################
#Table 4
#Total and Others Troop Fatalities
###################################



lmer_mdeaths<-lmer(ch_troops ~  lag_mdeaths*lag_lnpivotalaid*Contiguous + lag_n_contributors_troops + lag_mission_troops + lag_troops_mean + lag_dem + lag_lngdp + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_troops_mean>0, lag_mission_troops>0))
nobs(lmer_mdeaths)
summary(lmer_mdeaths)


lmer_odeaths<-lmer(ch_troops ~  lag_odeaths*lag_lnpivotalaid*Contiguous + lag_n_contributors_troops + lag_mission_troops + lag_troops_mean + lag_dem + lag_lngdp + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_troops_mean>0, lag_mission_troops>0))
nobs(lmer_odeaths)
summary(lmer_odeaths)


r.squaredGLMM(lmer_mdeaths)
r.squaredGLMM(lmer_odeaths)


#Additional marginal plots to visualize interactions 

# mdeaths_dat<-model.frame(lmer_mdeaths)
# cont_mdeaths_dat<-ggpredict(lmer_mdeaths, terms=c("lag_mdeaths", "Contiguous"))
# 
# plot(cont_mdeaths_dat, col="bw") + labs(x="Total Troop Deaths(t-1)", y="Change in Troops") + ggtitle("") + geom_rug(data=mdeaths_dat, aes(x=lag_mdeaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)
# 
# aid_mdeaths_dat<-ggpredict(lmer_mdeaths, terms=c("lag_mdeaths", "lag_lnpivotalaid[0, 10.77, 19.86]"))
# plot(aid_mdeaths_dat, col="bw") + labs(x="Total Troop Deaths(t-1)", y="Change in Troops") + ggtitle("") + geom_rug(data=mdeaths_dat, aes(x=lag_mdeaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)
# 
# aidcont_mdeaths_dat<-ggpredict(lmer_mdeaths, terms=c("lag_mdeaths", "lag_lnpivotalaid[0, 19.85]", "Contiguous"))
# plot(aidcont_mdeaths_dat, col="bw") + geom_line(aes(linetype=group)) + facet_wrap(~facet) + labs(x="Mission Deaths (t-1))", y="Change in Troops", colour="ln(Pivotal Aid(t-1))") + ggtitle("") + geom_rug(data=mdeaths_dat, aes(x=lag_mdeaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

# odeaths_dat<-model.frame(lmer_odeaths)
# cont_odeaths_dat<-ggpredict(lmer_odeaths, terms=c("lag_odeaths", "Contiguous"))
# 
# plot(cont_odeaths_dat, col="bw") + labs(x="Others' Troop Deaths(t-1)", y="Change in Troops") + ggtitle("") + geom_rug(data=odeaths_dat, aes(x=lag_odeaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)
# 
# aid_odeaths_dat<-ggpredict(lmer_odeaths, terms=c("lag_odeaths", "lag_lnpivotalaid[0, 10.77, 19.86]"))
# plot(aid_odeaths_dat, col="bw") + labs(x="Total Troop Deaths(t-1)", y="Change in Troops") + ggtitle("") + geom_rug(data=odeaths_dat, aes(x=lag_odeaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)
# 
# aidcont_odeaths_dat<-ggpredict(lmer_odeaths, terms=c("lag_odeaths", "lag_lnpivotalaid[0, 19.85]", "Contiguous"))
# plot(aidcont_odeaths_dat, col="bw") + geom_line(aes(linetype=group)) + facet_wrap(~facet) + labs(x="Others' Troop Deaths(t-1))", y="Change in Troops", colour="ln(Pivotal Aid(t-1))") + ggtitle("") + geom_rug(data=odeaths_dat, aes(x=lag_odeaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)



#Table 4 output
#Note variable labels, standard deviations for random effects, and r^2 manually inputted in ms manually 

stargazer(lmer_mdeaths, lmer_odeaths, title="Results from Multi-level Linear Regression", 
          align=T, digits=2, intercept.bottom=T, dep.var.labels="Change in Troops", star.cutoffs=c(0.05, 0.01, 0.001), no.space=T, style="io")



##################################
##########Appendix################
##################################

#########################
#Additional Descriptives
#########################

png('scatter_II.png', units="in", width=4.5, height=3.5, res=300)

ggplot(lmer_dat, aes(x=lag_troop_deaths, y=ch_troops)) + geom_point() + geom_smooth(method=lm) + labs(x="Troop Deaths(t-1)", y= "Change in Troops")

dev.off()

lmer_dat<- lmer_dat %>%
  mutate(lag_tfat_d=if_else(lag_troop_deaths>0,1,0))

lmer_dat<- lmer_dat %>% 
  mutate(withdraw=if_else(ch_troops<0, 1, 0)) %>%
  mutate(increase=if_else(ch_troops>0, 1, 0)) %>%
  mutate(nothing=if_else(ch_troops==0, 1, 0))

samp<- lmer_dat%>% 
  mutate(ch_type=case_when(
    ch_troops>0 ~ 1, 
    ch_troops==0 ~ 0, 
    ch_troops<0 ~ -1
  ))

samp<- samp %>%
  mutate(lag_tfat_d=if_else(lag_troop_deaths>0,1,0))

#Tabl 1A output
#Note variable labels, standard deviations for random effects, and r^2 manually inputted in ms manually 
withdraw_tab<- CrossTable(samp$lag_tfat_d, samp$ch_type, prop.c=FALSE, prop.t=FALSE, prop.chisq=FALSE)


inc<-samp %>%
  subset(ch_type==1)

dec<-samp %>%
  subset(ch_type==-1)

png('troopinc_II.png', units="in", width=4.5, height=3.5, res=300)

ggplot(inc, aes(x=ch_troops)) + geom_histogram(bins=15)  + xlab("Change in Troops")

dev.off()


png('troopdec_II.png', units="in", width=4.5, height=3.5, res=300)

ggplot(dec, aes(x=ch_troops)) + geom_histogram(bins=15) + xlab("Change in Troops")

dev.off()


##########################
#Checking Influential Obs
##########################

cooksd<-cooks.distance(lmer_full)

png('fig1a_II.png', units="in", width=4.5, height=3.5, res=300)

plot(cooksd, pch="*", cex=1, main="")
abline(h=0.0009638554, col="red")

dev.off()



samp$cooksd<-cooksd

influential<-as.numeric(names(cooksd)[(cooksd>0.0009638554)])
min(influential)

samp<-samp %>%
  mutate(outlier=ifelse(cooksd>(0.0009638554), 1, 0))

samp2<-samp%>%
  subset(outlier==0)

samp3<-samp%>%
  subset(outlier==1)


lm_f_outliers<-lmer(ch_troops ~ lag_troop_deaths*lag_lnpivotalaid*Contiguous + lag_dem +  lag_lngdp  + lag_troops_mean + lag_n_contributors_troops + lag_mission_troops + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(samp2, lag_troops_mean>0))
summary(lm_f_outliers)
nobs(lm_f_outliers)

#Table 2A output
#Note variable labels, standard deviations for random effects, and r^2 manually inputted in ms manually 

stargazer(lmer_full, lm_f_outliers, title="Results from Multi-level Linear Regression", 
          align=T, dep.var.labels="Change in Troops", no.space=T, star.cutoffs=0.05, style="io")

r.squaredGLMM(lmer_full)
r.squaredGLMM(lm_f_outliers)

summary(lmer_full)

#get predicted values

contdf_outliers<-ggpredict(lm_f_outliers, terms=c("lag_troop_deaths", "Contiguous"))
aiddf_outliers<-ggpredict(lm_f_outliers, terms=c("lag_troop_deaths", "lag_lnpivotalaid[1.69, 19.85]"))

##Contiguous Plot

png('fig2a_II.png', units="in", width=7, height=3.5, res=300)


plot(contdf_outliers, col="bw") + labs(x="Troop Deaths(t-1)", y="Change in Personnel") + scale_x_continuous(breaks=c(1,5,10), limits=c(0,10)) + ggtitle("") + geom_rug(data=samp2, aes(x=lag_troop_deaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()


##Pivotal Aid Plot


png('fig3a_II.png', units="in", width=7, height=3.5, res=300)

plot(aiddf_outliers, col="bw") + geom_line(aes(linetype=group)) + scale_x_continuous(breaks=c(1,5,10), limits=c(0,10)) + labs(x="Troop Deaths(t-1)", y="Change in Personnel", colour="ln(Pivotal Aid(t-1))") + ggtitle("") + geom_rug(data=samp2, aes(x=lag_troop_deaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()



############################################################

#Disaggregating Fatality Types

##Total deaths 

lm_total<-lmer(ch_total~  lag_total_deaths*lag_lnpivotalaid*Contiguous + lag_dem +  lag_lngdp  + lag_total_mean + lag_n_contributors_total + lag_mission_total + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_total_mean>0))
nobs(lm_total)
summary(lm_total)



all_dat<-model.frame(lm_total)


all_cont_df<-ggpredict(lm_total, terms=c("lag_total_deaths", "Contiguous"))
all_aid_df<-ggpredict(lm_total, terms=c("lag_total_deaths", "lag_lnpivotalaid[1.69, 19.85]"))


png('fig4A1_II.png', units="in", width=5, height=2.5, res=300)

plot(all_cont_df, col="bw") + labs(x="All PK Deaths(t-1)", y="Change in Personnel") + ggtitle("") + geom_rug(data=all_dat, aes(x=lag_total_deaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()

png('fig4A2_II.png', units="in", width=5, height=2.5, res=300)

plot(all_aid_df, col="bw") + labs(x="All PK Deaths(t-1)", y="Change in Personnel") + ggtitle("") + geom_rug(data=all_dat, aes(x=lag_total_deaths, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()


##Malicious attacks

lm_mal<-lmer(ch_total~  lag_Malicious*lag_lnpivotalaid*Contiguous + lag_dem +  lag_lngdp  + lag_total_mean + lag_n_contributors_total + lag_mission_total + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_total_mean>0))
summary(lm_mal)


mal_dat<-model.frame(lm_mal)


mal_cont_df<-ggpredict(lm_mal, terms=c("lag_Malicious", "Contiguous"))
mal_aid_df<-ggpredict(lm_mal, terms=c("lag_Malicious", "lag_lnpivotalaid[1.69, 19.85]"))

png('fig5A1_II.png', units="in", width=5, height=2.5, res=300)

plot(mal_cont_df, col="bw") + labs(x="Malicious Deaths(t-1)", y="Change in Personnel") + ggtitle("") + geom_rug(data=mal_dat, aes(x=lag_Malicious, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()

png('fig5A2_II.png', units="in", width=5, height=2.5, res=300)

plot(mal_aid_df, col="bw") + labs(x="Malicious Deaths(t-1)", y="Change in Personnel") + ggtitle("") + geom_rug(data=mal_dat, aes(x=lag_Malicious, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()


##Illness

lm_ill<-lmer(ch_total~  lag_Illness*lag_lnpivotalaid*Contiguous + lag_dem +  lag_lngdp  + lag_total_mean + lag_n_contributors_total + lag_mission_total + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_total_mean>0))
summary(lm_ill)


ill_dat<-model.frame(lm_ill)

ill_cont_df<-ggpredict(lm_ill, terms=c("lag_Illness", "Contiguous"))
ill_aid_df<-ggpredict(lm_ill, terms=c("lag_Illness", "lag_lnpivotalaid[1.69, 19.85]"))

png('fig6A1_II.png', units="in", width=5, height=2.5, res=300)

plot(ill_cont_df, col="bw") + labs(x="Illness Deaths(t-1)", y="Change in Personnel") + ggtitle("") + geom_rug(data=ill_dat, aes(x=lag_Illness, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()

png('fig6A2_II.png', units="in", width=5, height=2.5, res=300)

plot(ill_aid_df, col="bw") + labs(x="Illness Deaths(t-1)", y="Change in Personnel") + ggtitle("") + geom_rug(data=ill_dat, aes(x=lag_Illness, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()


##Accidents

lm_acc<-lmer(ch_total~  lag_Accident*lag_lnpivotalaid*Contiguous + lag_dem +  lag_lngdp  + lag_total_mean + lag_n_contributors_total + lag_mission_total + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_total_mean>0))
summary(lm_acc)



acc_dat<-model.frame(lm_acc)

acc_cont_df<-ggpredict(lm_acc, terms=c("lag_Accident", "Contiguous"))
acc_aid_df<-ggpredict(lm_acc, terms=c("lag_Accident", "lag_lnpivotalaid[1.69, 19.85]"))

png('fig7A1_II.png', units="in", width=5, height=2.5, res=300)

plot(acc_cont_df, col="bw") + labs(x="Accident Deaths(t-1)", y="Change in Personnel") + ggtitle("") + geom_rug(data=acc_dat, aes(x=lag_Accident, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()

png('fig7A2_II.png', units="in", width=5, height=2.5, res=300)

plot(acc_aid_df, col="bw") + labs(x="Accident Deaths(t-1)", y="Change in Personnel") + ggtitle("") + geom_rug(data=acc_dat, aes(x=lag_Accident, y=0), size=1, alpha=.25, position="jitter", sides="b", inherit.aes=F)

dev.off()

#Table 3A output
#Note variable labels, standard deviations for random effects, and r^2 manually inputted in ms manually 

stargazer(lm_total, lm_mal, lm_ill, lm_acc, title="Results from Multi-level Linear Regression", 
          align=T, dep.var.labels="Change in Personnel", no.space=T, star.cutoffs=0.05, style="io")


r.squaredGLMM(lm_total)
r.squaredGLMM(lm_mal)
r.squaredGLMM(lm_ill)
r.squaredGLMM(lm_acc)




#Other measures for 'intrinsic benefits'

#Refugees


lm_full_refugees<-lmer(ch_troops ~  lag_troop_deaths*lag_ln_refugees*lag_lnpivotalaid + lag_dem +  lag_lngdp  + lag_troops_mean + lag_n_contributors_troops + lag_mission_troops + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_troops_mean>0))
nobs(lm_full_refugees)
summary(lm_full_refugees)


png('fig8A_II.png', units="in", width=7, height=3.5, res=300)
refugees_plot<-plot_model(lm_full_refugees, type="pred", terms=c("lag_troop_deaths", "lag_lnpivotalaid[1.69, 19.85]", "lag_ln_refugees[0,14.53]"), title="", axis.title=c("Troop Deaths(t-1)", "Change in Troops"), mdrt.values="meansd", color="bw", se=lm_refugees_rse) + theme_bw()
refugees_plot
dev.off()


#Minimum Distance


lm_full_mindist<-lmer(ch_troops ~  lag_troop_deaths*ln_mindist*lag_lnpivotalaid + lag_n_contributors_troops + lag_mission_troops + lag_troops_mean + lag_dem + lag_lngdp + lag_ln_bdbest + (1|m_id) + (1|year) + (1|contributor_ccode), data=subset(APFP_dat, lag_troops_mean>0))
nobs(lm_full_mindist)
summary(lm_full_mindist)



png('fig9A_II.png', units="in", width=7, height=3.5, res=300)
mindist_plot<- plot_model(lm_full_mindist, type="pred", terms=c("lag_troop_deaths", "lag_lnpivotalaid[1.69, 19.85]", "ln_mindist[0,8.23,9.86]"), title="", axis.title=c("Troop Deaths(t-1)", "Change in Troops"), mdrt.values="meansd", color="bw", se=lm_mindist_rse) + theme_bw()
mindist_plot
dev.off()

aid_d_plot<- plot_model(lm_full_mindist, type="pred", terms=c("lag_troop_deaths", "lag_lnpivotalaid[1.69, 19.85]"), title="", axis.title=c("Troop Deaths(t-1)", "Change in Troops"), mdrt.values="meansd", color="bw", se=lm_mindist_rse) + theme_bw()
aid_d_plot

#Table 4A output
#Note variable labels, standard deviations for random effects, and r^2 manually inputted in ms manually 

stargazer(lm_full_refugees, lm_full_mindist, title="Results from Multi-level Linear Regression",
          omit.stat=c("adj.rsq", "f", "ser"), align=T, dep.var.labels="Change in Troops", star.cutoffs=0.05, no.space=T, star.char="*", style="io")



r.squaredGLMM(lm_full_refugees)
r.squaredGLMM(lm_full_mindist)

